# install.packages("readstata13")
rm(list=ls(all=TRUE))
library(readstata13)

# Set Working Directory
#setwd("")

dat<-read.dta13("FDI&Conflict.dta")

################################################################################
############### Calculate Marginal Effect: Model 2 in Table 1 ##################
################################################################################
mydata1<-dat%>%select(ccode, year, onset2v414, incidencev414, lwdist, cr_rfdicap,
                      lmrpe, llrgdpcap, llpop, llgrowthrate, lpolity2, 
                      lloil_gas_cap, elf85, relfrac, lmtnest, ncontig, coldwar, 
                      t1_onset2, t2_onset2, t3_onset2)%>%drop_na()

mydata1<-mydata1%>%mutate(s = as.numeric(incidencev414==1&onset2v414==0))%>%
  filter(s==0)

medianmat<-c(apply(mydata1[,6:20],2,median),1)

coef.sim<-read.dta13("sim_probit_m2.dta")
coefmat<-coef.sim[,1:16]

# Function for Calculating Marginal Effect
meffect<-function(medianmat,coefmat,N,x){
p1<- pnorm(t(as.matrix(medianmat))%*%t(coefmat))
p2<- pnorm(x*coefmat[,N]+t(as.matrix(medianmat))%*%t(coefmat))
print(c(mean(p2-p1),quantile(p2-p1,c(0.005,0.025,0.975,0.995))))
}

# Marginal Effect: FDI per Capita
meffect(medianmat,coefmat,1,sd(mydata1$cr_rfdicap))

################################################################################
#################   Figure 1: Plots of Marginal Effects ########################
################################################################################

# Model 1 in Table 4: RPE
coef.sim.rpe<-read.dta13("sim_rpe.dta")

capacity.sim<-seq(-1,2.5,0.01)
N<-length(capacity.sim)
p1<-matrix(NA,N,1000)
meffect<-matrix(NA, N, 5)
for (i in 1:N){
p1[i,]<-coef.sim.rpe$x1+coef.sim.rpe$x2*capacity.sim[i]
meffect[i,]<-quantile(p1[i,], c(0.025, 0.05, .5, 0.95, 0.975))
}


pdf(file="meffects_rpe.pdf",
    width=3.25,height=2.2)
par(mfrow=c(1,1),mar=c(2.1,2.2,.5,.8))
plot(capacity.sim, meffect[,3],ylim=c(-.02,0.07),type="l",ylab="Marginal Effects",
     col="black",lwd=.8 ,xlab="Relative Political Extraction (Mean Centered)",
     cex.lab=.6, cex.axis=.6,cex=.7,mgp = c(1.15,.5,0),axes=F)
box()
axis(1,at=seq(-1,2.5, by=0.5),cex=.6,cex.lab=.5,cex.axis=.6,mgp = c(1.0,.4,0))
axis(2,at=seq(-.02,.08,by=0.02),cex=.6,cex.lab=.5,cex.axis=.6,mgp = c(1.65,.5,0))
abline(h=0, lty=2,col="grey60",lwd=.8)
frame.y1 <- c(meffect[,1],rev(meffect[,5]))
frame.x1 <- c(capacity.sim, rev(capacity.sim))
polygon(frame.x1, frame.y1, density=-1, 
        col=rgb(166, 166, 166, 100, maxColorValue=255), border=NA)
dev.off()

#-------------------------------------------------------------------------------
# Model 2 in Table 4: Tax/GDP

coef.sim.tax<-read.dta13("sim_tax.dta")

capacity.sim<-seq(-2,1.66,0.01)
N<-length(capacity.sim)
p1<-matrix(NA,N,1000)
meffect<-matrix(NA, N, 5)
for (i in 1:N){
  p1[i,]<-coef.sim.tax$x1+coef.sim.tax$x2*capacity.sim[i]
  meffect[i,]<-quantile(p1[i,], c(0.025, 0.05, .5, 0.95, 0.975))
}


pdf(file="meffects_tax.pdf",
    width=3.25,height=2.2)
par(mfrow=c(1,1),mar=c(2.1,2.2,.5,.8))
plot(capacity.sim, meffect[,3],ylim=c(-.001,.09),type="l",ylab="Marginal Effects",
     col="black",lwd=.8, xlab="Tax (% of GDP, Mean Centered)",cex.lab=.6, 
     cex.axis=.6,cex=.7,mgp = c(1.15,.5,0),axes=F)
box()
axis(1,at=seq(-2,1.7, by=.5),cex=.6,cex.lab=.5,cex.axis=.6,mgp = c(1.0,.4,0))
axis(2,at=seq(0,.1,by=0.02),cex=.6,cex.lab=.5,cex.axis=.6,mgp = c(1.65,.5,0))
abline(h=0, lty=2,col="grey60",lwd=.8)
frame.y1 <- c(meffect[,1],rev(meffect[,5]))
frame.x1 <- c(capacity.sim, rev(capacity.sim))
polygon(frame.x1, frame.y1, density=-1, 
        col=rgb(166, 166, 166, 100, maxColorValue=255), border=NA)
dev.off()

#-------------------------------------------------------------------------------
# Model 3 in Table 4: Enrollment
coef.sim.pri<-read.dta13("sim_pri.dta")

capacity.sim<-seq(-.81, .73,0.01)
N<-length(capacity.sim)
p1<-matrix(NA,N,1000)
meffect<-matrix(NA, N, 5)
for (i in 1:N){
  p1[i,]<-coef.sim.pri$x1+coef.sim.pri$x2*capacity.sim[i]
  meffect[i,]<-quantile(p1[i,], c(0.025, 0.05, .5, 0.95, 0.975))
}

pdf(file="meffects_enroll.pdf",
    width=3.25,height=2.2)
par(mfrow=c(1,1),mar=c(2.1,2.2,.5,.8))
plot(capacity.sim, meffect[,3],ylim=c(-.01,0.081),type="l",ylab="Marginal Effects",
     col="black",lwd=.8, xlab="Primary School Enrollment Rate (Mean Centered)",
     cex.lab=.6, cex.axis=.6,cex=.7,mgp = c(1.15,.5,0),axes=F)
box()
axis(1,at=seq(-.8,.8, by=0.2),cex=.6,cex.lab=.5,cex.axis=.6,mgp = c(1.0,.4,0))
axis(2,at=seq(-.02,.08,by=0.02),cex=.6,cex.lab=.5,cex.axis=.6,mgp = c(1.65,.5,0))
abline(h=0, lty=2,col="grey60",lwd=.8)
frame.y1 <- c(meffect[,1],rev(meffect[,5]))
frame.x1 <- c(capacity.sim, rev(capacity.sim))
polygon(frame.x1, frame.y1, density=-1, 
        col=rgb(166, 166, 166, 100, maxColorValue=255), border=NA)
dev.off()
